#load datasets and libraries
library(rsample)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.4
## ✔ ggplot2   3.4.4     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(naniar)
library(corrplot)
## corrplot 0.92 loaded
Propdata <- read_csv("propofoldata.csv")
## Rows: 72 Columns: 44
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (23): Org_Name, Height, IBW, Triglyc_Max, CrCl, HR_Range_low, Min_Rass, ...
## dbl (21): Age, Sex, Weight, COVID, IVDA_Hx, Home_Opiod, Home_Benzo, Opiate_T...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(Propdata)
## # A tibble: 6 × 44
##   Org_Name     Age   Sex Height Weight IBW   COVID IVDA_Hx Home_Opiod Home_Benzo
##   <chr>      <dbl> <dbl> <chr>   <dbl> <chr> <dbl>   <dbl>      <dbl>      <dbl>
## 1 IU Health…    43     0 154.9   118.  47.7…     0       0          1          0
## 2 IU Health…    65     0 149.9    73.9 43.2…     0       0          0          0
## 3 IU Health…    64     0 180.3   128.  70.7…     1       0          0          0
## 4 IU Health…    67     0 157.5   155.  50.1…     0       0          0          0
## 5 IU Health…    63     0 167.6    79.4 59.2…     1       0          0          0
## 6 IU Health…    37     0 172.7   164   63.8…     0       0          1          1
## # ℹ 34 more variables: Opiate_Tol <dbl>, Triglyc_Max <chr>, Creatinine <dbl>,
## #   CrCl <chr>, Map_Range <dbl>, HR_Range_low <chr>, Total_days_prop <dbl>,
## #   Dur_prop_minutes <dbl>, Per_prop_over50 <dbl>, Max_Dose_Prop <dbl>,
## #   Min_Rass <chr>, Max_Rass <chr>, Dex_ord <chr>, Prop_dose_when_Dex <chr>,
## #   Dex_rate_when_prop <chr>, Prop_dec_afterDex <chr>, Fent_before_prop <chr>,
## #   Prop_dose_when_fent <chr>, Fent_max_when_prop <chr>,
## #   Prop_dec_after_fent <chr>, Numb_fent_bolus <chr>, Mode_CPOT <chr>, …
#show(Propdata)
#print(Propdata$Opiate_Tol)

Creating Analysis Dataset

IBW, COVID, IVDA_Hx, Triglyc_Max, Creatinine, CrCl, Map_Range, and HR_Range_low will all be dropped from the analysis dataset since they are not relevant to the analysis being completed.

#creating analysis dataset
df_eda <- select(Propdata, -c(IBW, COVID, IVDA_Hx,Triglyc_Max,Creatinine,CrCl,Map_Range,HR_Range_low))
df_eda[df_eda=="x"] <- NA
df_eda[df_eda=="X"] <- NA

x/X that represent NA/missing values is replaced with NA to be able to visualize and replace NA values.

#visualizing and appropriately replacing missing variables
vis_miss(df_eda[,1:12])

df_eda$Height <- ifelse(is.na(df_eda$Height), "Missing", df_eda$Height)
vis_miss(df_eda[,13:25])

df_eda$Min_Rass <- ifelse(is.na(df_eda$Min_Rass), "NotApp", df_eda$Min_Rass)
df_eda$Max_Rass <- ifelse(is.na(df_eda$Max_Rass), "NotApp", df_eda$Max_Rass)
df_eda$Dex_ord <- ifelse(is.na(df_eda$Dex_ord), "NotApp", df_eda$Dex_ord)
df_eda$Prop_dose_when_Dex <- ifelse(is.na(df_eda$Prop_dose_when_Dex), "NotApp", df_eda$Prop_dose_when_Dex)
df_eda$Dex_rate_when_prop <- ifelse(is.na(df_eda$Dex_rate_when_prop), "NotApp", df_eda$Dex_rate_when_prop)
df_eda$Prop_dec_afterDex <- ifelse(is.na(df_eda$Prop_dec_afterDex), "NotApp", df_eda$Prop_dec_afterDex)
df_eda$Fent_before_prop <- ifelse(is.na(df_eda$Fent_before_prop), "NotApp", df_eda$Fent_before_prop)
df_eda$Prop_dose_when_fent <- ifelse(is.na(df_eda$Prop_dose_when_fent), "NotApp", df_eda$Prop_dose_when_fent)
df_eda$Fent_max_when_prop <- ifelse(is.na(df_eda$Fent_max_when_prop), "NotApp", df_eda$Fent_max_when_prop)
df_eda$Prop_dec_after_fent <- ifelse(is.na(df_eda$Prop_dec_after_fent), "NotApp", df_eda$Prop_dec_after_fent)
df_eda$Numb_fent_bolus <- ifelse(is.na(df_eda$Numb_fent_bolus), "NotApp", df_eda$Numb_fent_bolus)
df_eda$Mode_CPOT <- ifelse(is.na(df_eda$Mode_CPOT), "NotApp", df_eda$Mode_CPOT)
df_eda$Was_ket_ord <- ifelse(is.na(df_eda$Was_ket_ord), "NotApp", df_eda$Was_ket_ord)
vis_miss(df_eda[,26:36])

df_eda$Prop_dose_when_ket <- ifelse(is.na(df_eda$Prop_dose_when_ket), "NotApp", df_eda$Prop_dose_when_ket)
df_eda$Ket_max_when_prop <- ifelse(is.na(df_eda$Ket_max_when_prop), "NotApp", df_eda$Ket_max_when_prop)
df_eda$Was_prop_dec <- ifelse(is.na(df_eda$Was_prop_dec), "NotApp", df_eda$Was_prop_dec)
df_eda$Con_Prop_NMBA <- ifelse(is.na(df_eda$Con_Prop_NMBA), "NotApp", df_eda$Con_Prop_NMBA)
vis_miss(df_eda)

View(df_eda)

Exploratory Data Analysis

Visualizing categorical variables

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=Org_Name))

Most patients are from the IU Health Ball Memorial Hospital. This is helpful since this is where we are concentrating our analysis.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Sex), fill=factor(Sex)))

There are more males than females.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Home_Opiod), fill=factor(Home_Opiod)))

Most patients do not have home opiod usage.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Home_Benzo), fill=factor(Home_Benzo)))

Most patients do not have home benzo usage.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Opiate_Tol), fill=factor(Opiate_Tol)))

Most patients have no tolerance for Opiate.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Min_Rass), fill=factor(Min_Rass)))

Min RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients minimum RASS when propofol rate >50 mcg/kg/min was that they were unarousable. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Max_Rass), fill=factor(Max_Rass)))

Max RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients maximum RASS when propofol rate >50 mcg/kg/min was that they would pull at or attempt to remove tubes alongside generally aggressive behavior. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Dex_ord), fill=factor(Dex_ord)))

Dexmedetomidine was ordered before propofol rate reached >50 mcg/kg/min at any point. Most patients were not ordered dexmedetomidine.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Prop_dec_afterDex), fill=factor(Prop_dec_afterDex)))

Propofol was decreased to <50 mcg/kg/min after dexmedetomidine was added (did not go back above 50). Since most patients were not ordered dexmedetomidine, this is not applicable to most patients.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Fent_before_prop), fill=factor(Fent_before_prop)))

Fentanyl was ordered before propofol rate was >50 mcg/kg/min. Most patients did have fentanyl ordered for them before the propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Prop_dec_after_fent), fill=factor(Prop_dec_after_fent)))

Propofol was decreased to <50 mcg/kg/min after fentanyl was added. This was not applicable to most patients.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Mode_CPOT), fill=factor(Mode_CPOT)))

Mode CPOT when propofol rate >50 mcg/kg/min . “The CPOT is a behavioural assessment pain scale for patients unable to verbalise pain.” (2) Most patients had a mode CPOT of 0, which means that overall, the patient was relaxed, had no muscle tension, and had an absence of movement. Link to learn more about the scale: https://ccs-sth.org/resources/Documents/Sedation%20Analgesia%20Delirium%20in%20CC/CCS%20STH%20Guideline%20template%20CPOT.pdf

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Was_ket_ord), fill=factor(Was_ket_ord)))

Most patients did not have ketamine prescribed to them before propofol rate was >50 mcg/kg/min and out of the ones who did, most did not have it ordered.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Was_prop_dec), fill=factor(Was_prop_dec)))

Since most patients did not have ketamine prescribed to them, it was not applicable to most patients if propofol decreased to <50 mcg/kg/min after ketamine was added or the dose was increased.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Rec_pain_med), fill=factor(Rec_pain_med)))

Most patients did not receive any scheduled pain medications when propofol >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(On_benzo), fill=factor(On_benzo)))

Most patients were not on scheduled benzo.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Vas), fill=factor(Con_Prop_Vas)))

More patients did not receive concurrent propofol and vasopressor when propofol >50 mcg/kg/min than those who did.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Dex), fill=factor(Con_Prop_Dex)))

More patients did not receive concurrent propofol and dexmedetomidine than those who did.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Fen), fill=factor(Con_Prop_Fen)))

More patients did receive concurrent propofol and fentanyl infusion than those who did not.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Ket), fill=factor(Con_Prop_Ket)))

Most patients did not receive concurrent propofol and ketamine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Ben), fill=factor(Con_Prop_Ben)))

Most patients did not receive concurrent propofol and benzodiazepine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_NMBA), fill=factor(Con_Prop_NMBA)))

More patients did not receive concurrent propofol and NMBA infusion than those who did.

Visualizing covariation between two categorical variables

Creating initialisms of Org_Name values to make bar graphs easier to interpret.

df_eda$Org_Name[df_eda$Org_Name=="IU Health Arnett Hospital"] <- 'IU HAH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Ball Memorial Hospital"] <- 'IU HBMH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Bloomington Hospital"] <- 'IU HBloomH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Methodist Hospital"] <- 'IU HMH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health University Hospital"] <- 'IU HUH'

Initialisms for hospital names:

  • IU Health Arnett Hospital <- ‘IU HAH’
  • IU Health Ball Memorial Hospital <- ‘IU HBMH’
  • IU Health Bloomington Hospital <- ‘IU HBloomH’
  • IU Health Methodist Hospital <- ‘IU HMH’
  • IU Health University Hospital <- ‘IU HUH’
ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Opiod)), position = "dodge")

The organization with the most patients that have home opiod usage is IU HBMH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Benzo)), position = "dodge")

The organization that has the most patients with home benzo usage is IU HBloomH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Opiate_Tol)), position = "dodge")

IU HBMH has the most patients with high tolerance of opiates.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_afterDex)), position = "dodge")

Only IU HBMH and IU HBloomH had propofol decreases to <50 mcg/kg/min after dexmedetomidine was added.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Fent_before_prop)), position = "dodge")

IU HBMH had the highest number of fentanyl orders before propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_after_fent)), position = "dodge")

Only IU HBloomH and IU HMH had propofol decreases to <50 mcg/kg/min after fentanyl was added.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_ket_ord)), position = "dodge")

IU HMH had the highest number of ketamine orders before propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_prop_dec)), position = "dodge")

Propofol was only decreased to <50 mcg/kg/min after ketamine was added or dose increased at IU HBMH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Rec_pain_med)), position = "dodge")

Only patients at IU HMH and IU HUH received scheduled pain medications when propofol >50.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Vas)), position = "dodge")

Patients received concurrent propofol and vasopressor when propofol >50 mcg/kg/min at IU HBMH, IU HMH, and IU HBloomH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Dex)), position = "dodge")

Only IU HUH did not have patients who received concurrent propofol and dexmedetomidine.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Fen)), position = "dodge")

Concurrent propofol and fentanyl infusion is a popular choice at all organizations except IU HUH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ket)), position = "dodge")

On the contrary, concurrent propofol and ketamine infusion was not a popular choice at the organizations.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ben)), position = "dodge")

IU HMH had the most patients that received concurrent propofol and benzodiazepine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_NMBA)), position = "dodge")

IU HBMH had the most patients that received concurrent propofol and NMBA infusion.

Visualizing continuous variables

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Age), binwidth = 3)

Most patients are between 50 and 70 years old.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Height)), binwidth = 10)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Most patients are around 180cm in height (roughly 5 ft 11in).

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Weight), binwidth = 5)

Weights are pretty evenly distributed with a few obvious outliers.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Total_days_prop), binwidth = 1)

Most patients were on propofol for 15 days or less, with a few outliers towards 20 days and 30 days.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Dur_prop_minutes), binwidth = 2000)

Most patients were on propofol for less than 10000 minutes (which is roughly less than a week) with a few higher outliers.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Per_prop_over50), binwidth = 5)

The % of time of propofol above 50 mcg/kg/min is relatively spread across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Max_Dose_Prop), binwidth = 3)

The max dose of propofol in mcg/kg/min greatly varies across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_Dex)), binwidth = 3)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).

This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when dexmedetomidine was added varied across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Dex_rate_when_prop)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the dexmedetomidine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was usually around 0.5 mcg/kg/hr.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_fent)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 59 rows containing non-finite values (`stat_bin()`).

This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was varied.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Fent_max_when_prop)), binwidth = 5)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 22 rows containing non-finite values (`stat_bin()`).

The fentanyl max rate in mcg/hr when propofol >50 mcg/kg/min was between 100-150 mcg/hr for most patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Numb_fent_bolus)), binwidth = 5)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 13 rows containing non-finite values (`stat_bin()`).

Most patients did not receive fentanyl bolus when propofol rate >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_ket)), binwidth = 3)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 65 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was usually 50 mcg/kg/min or above.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Ket_max_when_prop)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 61 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the ketamine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was spread between 0.0 to 2.0 mcg/kg/hr.

Visualizing covariation between a categorical variable and a continuous variable

The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. Based off evaluating the EDA results, the selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA. Below are some additional visualizations showing covariation between per_prop_over50 and the selected covariates.

ggplot(data = df_eda, mapping = aes(x = factor(Org_Name), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Home_Opiod), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Home_Benzo), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Opiate_Tol), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Min_Rass), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Max_Rass), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Dex_ord), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Prop_dec_afterDex), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Fent_before_prop), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Mode_CPOT), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Was_ket_ord), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Was_prop_dec), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Rec_pain_med), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(On_benzo), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Vas), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Dex), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Fen), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ket), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ben), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_NMBA), y = Per_prop_over50)) + geom_boxplot()

Fitting GLM by selecting a quantitative response variable column and a set of covariates

The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. The selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA.

glm_mod1 <- glm(Per_prop_over50~ Home_Opiod+Opiate_Tol+Dex_ord+Prop_dose_when_Dex+Fent_before_prop+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+On_benzo+Con_Prop_Vas+Con_Prop_Dex+Con_Prop_Ket+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod1)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Home_Opiod + Opiate_Tol + Dex_ord + 
##     Prop_dose_when_Dex + Fent_before_prop + Prop_dec_after_fent + 
##     Mode_CPOT + Was_ket_ord + Rec_pain_med + On_benzo + Con_Prop_Vas + 
##     Con_Prop_Dex + Con_Prop_Ket + Con_Prop_Ben + Con_Prop_NMBA, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (3 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                11.0772    26.0465   0.425  0.67366   
## Home_Opiod                 11.4687    10.0100   1.146  0.26097   
## Opiate_Tol                  3.7031     6.8949   0.537  0.59517   
## Dex_ord1                   17.4662    41.1774   0.424  0.67447   
## Dex_ordNotApp              25.4739    40.2750   0.632  0.53185   
## Prop_dose_when_Dex10       39.3687    34.3706   1.145  0.26109   
## Prop_dose_when_Dex100      29.1390    23.9251   1.218  0.23275   
## Prop_dose_when_Dex21.2404  28.3292    28.2648   1.002  0.32423   
## Prop_dose_when_Dex40      -11.1913    23.8410  -0.469  0.64217   
## Prop_dose_when_Dex50       22.6597    24.5651   0.922  0.36366   
## Prop_dose_when_Dex50.0002   7.8072    23.2330   0.336  0.73918   
## Prop_dose_when_Dex51.2821  61.9114    57.6052   1.075  0.29105   
## Prop_dose_when_Dex52       31.9668    36.3055   0.880  0.38559   
## Prop_dose_when_Dex53.8642  45.4851    30.6342   1.485  0.14803   
## Prop_dose_when_Dex54.1463  32.0240    36.1578   0.886  0.38284   
## Prop_dose_when_Dex59.991   27.6077    24.5027   1.127  0.26879   
## Prop_dose_when_Dex60       21.7064    29.6990   0.731  0.47052   
## Prop_dose_when_Dex65        9.7494    23.9251   0.407  0.68654   
## Prop_dose_when_Dex70        2.7771    21.2857   0.130  0.89707   
## Prop_dose_when_Dex80       11.2612    34.7263   0.324  0.74797   
## Prop_dose_when_Dex89.8261  26.3817    25.4026   1.039  0.30732   
## Prop_dose_when_DexNotApp  -14.8992    43.0625  -0.346  0.73177   
## Fent_before_prop1          -5.6738    14.8100  -0.383  0.70434   
## Fent_before_propNotApp     10.1541    15.8590   0.640  0.52686   
## Prop_dec_after_fent1      -19.2199    17.0755  -1.126  0.26927   
## Prop_dec_after_fentNotApp   8.4598    16.7658   0.505  0.61753   
## Mode_CPOT1                  3.4145    24.2573   0.141  0.88900   
## Mode_CPOT2                 11.4906    12.1827   0.943  0.35312   
## Mode_CPOT3                 14.5515    11.0369   1.318  0.19733   
## Mode_CPOT4                  0.8971    12.7054   0.071  0.94418   
## Mode_CPOT5                 24.4714     9.8122   2.494  0.01837 * 
## Mode_CPOT6                 37.9780    13.7865   2.755  0.00989 **
## Mode_CPOT7                      NA         NA      NA       NA   
## Mode_CPOT8                -14.0736    26.9472  -0.522  0.60532   
## Mode_CPOTNotApp            -3.6458    17.8405  -0.204  0.83945   
## Was_ket_ord1              -12.6255    23.6656  -0.533  0.59762   
## Was_ket_ordNotApp         -23.6351    19.8455  -1.191  0.24301   
## Rec_pain_med              -14.5930    13.7012  -1.065  0.29533   
## On_benzo                   -8.1254    18.2073  -0.446  0.65861   
## Con_Prop_Vas                4.5572     8.1382   0.560  0.57965   
## Con_Prop_Dex                    NA         NA      NA       NA   
## Con_Prop_Ket                    NA         NA      NA       NA   
## Con_Prop_Ben              -17.5800     9.6445  -1.823  0.07831 . 
## Con_Prop_NMBA1             10.8414     6.6535   1.629  0.11368   
## Con_Prop_NMBANotApp         8.5955    24.0118   0.358  0.72287   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 287.5668)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance:  8627  on 30  degrees of freedom
## AIC: 634.92
## 
## Number of Fisher Scoring iterations: 2

Creating new model with the removal of home_opiod, opiate_tol, dex_ord, fent_before_prop, on_benzo, con_prop_vas, con_prop_dex, and con_prop_ket due to insignificance.

glm_mod2 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod2)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT + Was_ket_ord + Rec_pain_med + Con_Prop_Ben + Con_Prop_NMBA, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                  3.464     19.087   0.181  0.85694   
## Prop_dose_when_Dex10        39.541     31.496   1.255  0.21700   
## Prop_dose_when_Dex100       30.581     22.727   1.346  0.18641   
## Prop_dose_when_Dex21.2404   39.896     23.530   1.696  0.09815 . 
## Prop_dose_when_Dex40        10.636     22.727   0.468  0.64247   
## Prop_dose_when_Dex50         9.453     22.667   0.417  0.67900   
## Prop_dose_when_Dex50.0002   16.457     22.667   0.726  0.47228   
## Prop_dose_when_Dex51.2821   63.853     36.876   1.732  0.09146 . 
## Prop_dose_when_Dex52        54.704     32.230   1.697  0.09782 . 
## Prop_dose_when_Dex53.8642   46.008     25.953   1.773  0.08429 . 
## Prop_dose_when_Dex54.1463   59.457     27.126   2.192  0.03459 * 
## Prop_dose_when_Dex59.991    34.877     22.157   1.574  0.12376   
## Prop_dose_when_Dex60        39.255     24.206   1.622  0.11314   
## Prop_dose_when_Dex65        11.191     22.727   0.492  0.62525   
## Prop_dose_when_Dex70        18.863     19.250   0.980  0.33334   
## Prop_dose_when_Dex80        32.804     30.175   1.087  0.28382   
## Prop_dose_when_Dex89.8261   28.054     22.157   1.266  0.21317   
## Prop_dose_when_DexNotApp    20.246     14.433   1.403  0.16880   
## Prop_dec_after_fent1       -21.775     15.536  -1.402  0.16917   
## Prop_dec_after_fentNotApp   -1.862     10.634  -0.175  0.86194   
## Mode_CPOT1                  -4.981     14.812  -0.336  0.73849   
## Mode_CPOT2                   3.884     11.433   0.340  0.73595   
## Mode_CPOT3                   9.856      9.335   1.056  0.29774   
## Mode_CPOT4                  10.898     10.277   1.060  0.29565   
## Mode_CPOT5                  28.093      8.862   3.170  0.00301 **
## Mode_CPOT6                  32.399     13.153   2.463  0.01841 * 
## Mode_CPOT7                      NA         NA      NA       NA   
## Mode_CPOT8                  -4.369     23.829  -0.183  0.85549   
## Mode_CPOTNotApp            -11.875     17.618  -0.674  0.50436   
## Was_ket_ord1                 5.466     19.175   0.285  0.77714   
## Was_ket_ordNotApp           -9.113     14.031  -0.650  0.51992   
## Rec_pain_med                -6.995     12.754  -0.548  0.58661   
## Con_Prop_Ben               -13.586      8.041  -1.690  0.09931 . 
## Con_Prop_NMBA1               8.191      5.391   1.520  0.13691   
## Con_Prop_NMBANotApp          6.858     22.667   0.303  0.76389   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 295.4365)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 11227  on 38  degrees of freedom
## AIC: 637.88
## 
## Number of Fisher Scoring iterations: 2

Creating a new model with the removal of was_ket_ord and rec_pain_med due to insignificance.

glm_mod3 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod3)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT + Con_Prop_Ben + Con_Prop_NMBA, family = "gaussian", 
##     data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                 -5.200     14.034  -0.371  0.71290   
## Prop_dose_when_Dex10        38.119     30.742   1.240  0.22204   
## Prop_dose_when_Dex100       31.309     22.416   1.397  0.17001   
## Prop_dose_when_Dex21.2404   37.383     23.138   1.616  0.11384   
## Prop_dose_when_Dex40        11.364     22.416   0.507  0.61491   
## Prop_dose_when_Dex50         9.381     22.387   0.419  0.67739   
## Prop_dose_when_Dex50.0002   16.529     22.387   0.738  0.46453   
## Prop_dose_when_Dex51.2821   62.910     27.929   2.253  0.02971 * 
## Prop_dose_when_Dex52        48.761     28.948   1.684  0.09969 . 
## Prop_dose_when_Dex53.8642   60.659     22.387   2.710  0.00979 **
## Prop_dose_when_Dex54.1463   69.298     22.416   3.091  0.00357 **
## Prop_dose_when_Dex59.991    35.981     21.834   1.648  0.10700   
## Prop_dose_when_Dex60        39.047     23.901   1.634  0.10998   
## Prop_dose_when_Dex65        11.919     22.416   0.532  0.59778   
## Prop_dose_when_Dex70        22.763     18.144   1.255  0.21675   
## Prop_dose_when_Dex80        43.158     26.002   1.660  0.10459   
## Prop_dose_when_Dex89.8261   29.158     21.834   1.335  0.18908   
## Prop_dose_when_DexNotApp    20.831     14.252   1.462  0.15146   
## Prop_dec_after_fent1       -20.944     15.054  -1.391  0.17165   
## Prop_dec_after_fentNotApp   -3.039     10.281  -0.296  0.76906   
## Mode_CPOT1                  -1.902     12.743  -0.149  0.88206   
## Mode_CPOT2                   1.866     11.114   0.168  0.86747   
## Mode_CPOT3                   9.024      9.122   0.989  0.32830   
## Mode_CPOT4                  10.521     10.082   1.044  0.30281   
## Mode_CPOT5                  29.029      8.541   3.399  0.00152 **
## Mode_CPOT6                  31.886     12.971   2.458  0.01828 * 
## Mode_CPOT7                      NA         NA      NA       NA   
## Mode_CPOT8                   2.301     18.826   0.122  0.90331   
## Mode_CPOTNotApp            -11.732     17.354  -0.676  0.50279   
## Con_Prop_Ben               -11.000      7.698  -1.429  0.16060   
## Con_Prop_NMBA1               8.847      5.268   1.679  0.10068   
## Con_Prop_NMBANotApp          6.785     22.387   0.303  0.76336   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 288.4435)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 11826  on 41  degrees of freedom
## AIC: 635.63
## 
## Number of Fisher Scoring iterations: 2

Creating a new model with the removal of con_prop_ben and con_prop_NMBA due to insignificance.

glm_mod4 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT ,data=df_eda,family="gaussian")
summary(glm_mod4)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT, family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                 0.6689    12.0325   0.056  0.95592   
## Prop_dose_when_Dex10       35.6161    30.6364   1.163  0.25128   
## Prop_dose_when_Dex100      27.5037    22.1297   1.243  0.22051   
## Prop_dose_when_Dex21.2404  31.4240    22.1297   1.420  0.16266   
## Prop_dose_when_Dex40        7.5586    22.1297   0.342  0.73431   
## Prop_dose_when_Dex50        3.4673    20.4161   0.170  0.86592   
## Prop_dose_when_Dex50.0002  21.5706    22.1297   0.975  0.33502   
## Prop_dose_when_Dex51.2821  60.8641    27.8911   2.182  0.03448 * 
## Prop_dose_when_Dex52       57.7160    27.8911   2.069  0.04442 * 
## Prop_dose_when_Dex53.8642  65.7009    22.1297   2.969  0.00482 **
## Prop_dose_when_Dex54.1463  65.4926    22.1297   2.959  0.00495 **
## Prop_dose_when_Dex59.991   32.1317    20.9987   1.530  0.13313   
## Prop_dose_when_Dex60       36.0523    23.7780   1.516  0.13662   
## Prop_dose_when_Dex65        8.1141    22.1297   0.367  0.71563   
## Prop_dose_when_Dex70       22.7545    17.5799   1.294  0.20230   
## Prop_dose_when_Dex80       32.2657    25.1141   1.285  0.20560   
## Prop_dose_when_Dex89.8261  25.3087    20.9987   1.205  0.23455   
## Prop_dose_when_DexNotApp   18.7854    13.7528   1.366  0.17890   
## Prop_dec_after_fent1      -22.5511    14.5088  -1.554  0.12727   
## Prop_dec_after_fentNotApp  -5.1024     9.9406  -0.513  0.61031   
## Mode_CPOT1                 -4.7385    12.5270  -0.378  0.70706   
## Mode_CPOT2                 -1.3287    10.3857  -0.128  0.89878   
## Mode_CPOT3                  8.9382     9.1286   0.979  0.33286   
## Mode_CPOT4                 10.5658     9.6682   1.093  0.28041   
## Mode_CPOT5                 28.2192     8.5576   3.298  0.00194 **
## Mode_CPOT6                 38.9734    12.5270   3.111  0.00327 **
## Mode_CPOT7                      NA         NA      NA       NA   
## Mode_CPOT8                -10.4584    17.4390  -0.600  0.55178   
## Mode_CPOTNotApp           -13.4918    17.4390  -0.774  0.44327   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 294.3858)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 12953  on 44  degrees of freedom
## AIC: 636.18
## 
## Number of Fisher Scoring iterations: 2

Performing variable selection using deviance tables

anova(glm_mod4,test="Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Per_prop_over50
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                   71      31735            
## Prop_dose_when_Dex  17  10969.2        54      20766 0.003101 **
## Prop_dec_after_fent  2   1061.9        52      19704 0.164695   
## Mode_CPOT            8   6750.8        44      12953 0.003453 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The final model for the % of time of propofol above 50 mcg/kg/min (Per_prop_over50) based on the deviance table can be written as follows:

Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT

glm_mod5 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Mode_CPOT,data=df_eda,family="gaussian")
summary(glm_mod5)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Mode_CPOT, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                 1.9746    11.8023   0.167  0.86786   
## Prop_dose_when_Dex10       10.8880    27.1926   0.400  0.69071   
## Prop_dose_when_Dex100      21.0956    20.9709   1.006  0.31970   
## Prop_dose_when_Dex21.2404  25.0159    20.9709   1.193  0.23903   
## Prop_dose_when_Dex40        1.1504    20.9709   0.055  0.95649   
## Prop_dose_when_Dex50        6.6714    20.2591   0.329  0.74342   
## Prop_dose_when_Dex50.0002  15.1625    20.9709   0.723  0.47332   
## Prop_dose_when_Dex51.2821  53.5848    27.1926   1.971  0.05481 . 
## Prop_dose_when_Dex52       50.4366    27.1926   1.855  0.07004 . 
## Prop_dose_when_Dex53.8642  59.2928    20.9709   2.827  0.00693 **
## Prop_dose_when_Dex54.1463  59.0844    20.9709   2.817  0.00711 **
## Prop_dose_when_Dex59.991   30.2333    20.2591   1.492  0.14244   
## Prop_dose_when_Dex60       28.3028    22.4704   1.260  0.21418   
## Prop_dose_when_Dex65        1.7060    20.9709   0.081  0.93552   
## Prop_dose_when_Dex70       18.2269    17.4245   1.046  0.30100   
## Prop_dose_when_Dex80       24.9863    24.2733   1.029  0.30869   
## Prop_dose_when_Dex89.8261  23.4104    20.2591   1.156  0.25383   
## Prop_dose_when_DexNotApp   11.5061    11.7672   0.978  0.33328   
## Mode_CPOT1                 -3.8672    12.6234  -0.306  0.76072   
## Mode_CPOT2                 -0.4575    10.4532  -0.044  0.96528   
## Mode_CPOT3                  9.8094     9.1777   1.069  0.29072   
## Mode_CPOT4                  6.0560     9.3834   0.645  0.52188   
## Mode_CPOT5                 29.5605     8.6006   3.437  0.00126 **
## Mode_CPOT6                 39.8446    12.6234   3.156  0.00282 **
## Mode_CPOT7                      NA         NA      NA       NA   
## Mode_CPOT8                 -9.5872    17.5953  -0.545  0.58847   
## Mode_CPOTNotApp           -12.6206    17.5953  -0.717  0.47683   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 300.4846)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 13822  on 46  degrees of freedom
## AIC: 636.86
## 
## Number of Fisher Scoring iterations: 2

Running diagnostics and checking model adequacy

par(mfrow=c(2,2))
plot(glm_mod5)
## Warning: not plotting observations with leverage one:
##   1, 2, 4, 7, 8, 9, 12, 17, 18, 19, 20, 21, 22, 26, 71
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Evaluating evidence of interaction between a quantitative and a qualitative variable

There is evidence of interaction between a quantitative predictor and a qualitative predictor based off the EDA. There seems to be some evidence of interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative). An ANCOVA will be run.

#Starting with ANOVA
Propdata_anova = select(df_eda, Per_prop_over50,Org_Name, Total_days_prop)
glm_anova <- glm(Per_prop_over50~ factor(Org_Name),data=Propdata_anova,family="gaussian")
summary(glm_anova)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name), family = "gaussian", 
##     data = Propdata_anova)
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                  11.047      9.426   1.172    0.245  
## factor(Org_Name)IU HBloomH    3.073     11.329   0.271    0.787  
## factor(Org_Name)IU HBMH      19.147      9.886   1.937    0.057 .
## factor(Org_Name)IU HMH       -4.314     10.539  -0.409    0.684  
## factor(Org_Name)IU HUH       -2.468     14.399  -0.171    0.864  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 355.4004)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 23812  on 67  degrees of freedom
## AIC: 634.02
## 
## Number of Fisher Scoring iterations: 2

There is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.

#Running ANCOVA
glm_ancova <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop, data=Propdata_anova,family="gaussian")
summary(glm_ancova)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop, 
##     family = "gaussian", data = Propdata_anova)
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                11.64986   10.11925   1.151   0.2538  
## factor(Org_Name)IU HBloomH  2.98745   11.42234   0.262   0.7945  
## factor(Org_Name)IU HBMH    19.00153    9.99409   1.901   0.0616 .
## factor(Org_Name)IU HMH     -4.36904   10.62054  -0.411   0.6821  
## factor(Org_Name)IU HUH     -2.65121   14.54262  -0.182   0.8559  
## Total_days_prop            -0.05482    0.31810  -0.172   0.8637  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 360.623)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 23801  on 66  degrees of freedom
## AIC: 635.99
## 
## Number of Fisher Scoring iterations: 2
anova(glm_ancova)
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Per_prop_over50
## 
## Terms added sequentially (first to last)
## 
## 
##                  Df Deviance Resid. Df Resid. Dev
## NULL                                71      31735
## factor(Org_Name)  4   7923.1        67      23812
## Total_days_prop   1     10.7        66      23801

Once again, there is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.

#interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative)
par(mar=c(1,1,1,1))
par(mfrow=c(1,1))
ggplot(data = Propdata_anova, aes(x=Total_days_prop, y=Per_prop_over50, color=factor(Org_Name)))+
  geom_point() + geom_smooth() + xlab("total days on propofol infusion (Total_days_prop)") + ylab("% of time of propofol above 50 mcg/kg/min (Per_prop_over50)")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small.  fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4.93
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 8.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 145.68
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.  fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 4.93
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 8.07
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 145.68
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1.1299e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 121
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 4
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 1.1299e-17
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 121
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small.  fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 3.97
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 5.03
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1.0609
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.  fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 3.97
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 5.03
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1.0609
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

#testing for interaction
glm_ancova1 <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop +factor(Org_Name)*Total_days_prop, data=Propdata_anova, family="gaussian")
summary(glm_ancova1)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop + 
##     factor(Org_Name) * Total_days_prop, family = "gaussian", 
##     data = Propdata_anova)
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 20.9435    21.5036   0.974    0.334
## factor(Org_Name)IU HBloomH                   2.7782    23.5070   0.118    0.906
## factor(Org_Name)IU HBMH                      6.6013    21.9946   0.300    0.765
## factor(Org_Name)IU HMH                     -12.6945    23.2714  -0.545    0.587
## factor(Org_Name)IU HUH                     -19.5092    40.3735  -0.483    0.631
## Total_days_prop                             -0.8997     1.7499  -0.514    0.609
## factor(Org_Name)IU HBloomH:Total_days_prop  -0.1170     1.9014  -0.062    0.951
## factor(Org_Name)IU HBMH:Total_days_prop      1.2169     1.7991   0.676    0.501
## factor(Org_Name)IU HMH:Total_days_prop       0.7481     1.9037   0.393    0.696
## factor(Org_Name)IU HUH:Total_days_prop       1.8315     4.5654   0.401    0.690
## 
## (Dispersion parameter for gaussian family taken to be 367.4723)
## 
##     Null deviance: 31735  on 71  degrees of freedom
## Residual deviance: 22783  on 62  degrees of freedom
## AIC: 640.84
## 
## Number of Fisher Scoring iterations: 2
#running model diagnostics
anova(glm_ancova, glm_ancova1)
## Analysis of Deviance Table
## 
## Model 1: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop
## Model 2: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop + factor(Org_Name) * 
##     Total_days_prop
##   Resid. Df Resid. Dev Df Deviance
## 1        66      23801            
## 2        62      22783  4   1017.8

The interaction term is significant as the difference in deviance between the larger model and the smaller model is statistically significant.

Final model with interaction term is as follows:

Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT - 0.518 * Total_days_prop

The final explanatory variables for the % of time of propofol above 50 mcg/kg/min are Prop_dose_when_dex, Mode_CPOT, and Total_days_prop.

Bibliography

  1. Physiopedia. (n.d.). Richmond agitation-sedation scale (RASS). Physiopedia. https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS) Description of Richmond Agitation-Sedation Scale (RASS).

  2. South Tees Hospitals NHS FT. (2023, January 20). CPOT = Critical-Care Pain Observation Tool. Critical Care Services - The James Cook University Hospital. https://ccs-sth.org/resources/Documents/Sedation Analgesia Delirium in CC/CCS STH Guideline template CPOT.pdf Description of CPOT = Critical-Care Pain Observation Tool.